% fig2b2.m ——— 并行化绘制 LLE vs ω（调用 LEs_analytic）
clear; clc; close all;

%% —— 固定模型参数 —— 
a     = 0.5;   b     = 0.215;
beta  = 0.2;   A     = 0.95;
alpha = 1.8;

%% —— 仿真时间设置 —— 
dt = 0.01;  T = 200;

%% —— 扫描 ω 列表 —— 
omega_list = linspace(0.4, 0.8, 1000);
N = numel(omega_list);
lle_vals = zeros(1, N);
X0 = [0.01; 0.02; 0.02];

%% —— 并行计算 LLE —— 
parfor k = 1:N
    omega_k = omega_list(k);
    f   = @(t,X) rhs(t,   X, a, b, beta, A, omega_k, alpha);
    j = @(X)   jac(X, a, b, beta, A, omega_k, alpha);
    lle_vals(k) = LEs(f, j, X0, dt, T);
end

%% —— 绘图 —— 
figure;
plot(omega_list, lle_vals, '.-','MarkerSize',6); hold on;
yline(0,'r--','LineWidth',1.2);
xlabel('\omega','FontSize',12);
ylabel('LE 指数','FontSize',12);
title('LE 随 \omega 变化 (\alpha = 1.8)','FontSize',14);
grid on;

% % fig2b2_noise.m ——— 加入白噪声的 LLE vs ω 曲线
% clear; clc; close all;
% 
% %% —— 模型参数 —— 
% a     = 0.5;
% b     = 0.215;
% beta  = 0.2;
% A     = 0.95;
% alpha = 1.8;         % 固定 α
% mu_s  = 0;           % μ_s 未指定则设为 0
% 
% %% —— 噪声参数 —— 
% D     = 0.01;        % 噪声强度
% dt    = 0.01;
% T     = 200;
% warmup = 100;
% 
% omega_list = linspace(0.4, 0.8, 200);
% N = numel(omega_list);
% lle_vals = zeros(1,N);
% 
% X0 = [0.01; 0.02; 0.02];
% 
% %% —— 启动并行池 —— 
% if isempty(gcp('nocreate')), parpool; end
% 
% %% —— 并行计算 —— 
% parfor k = 1:N
%     omega_k = omega_list(k);
% 
%     % —— 构造噪声项 Sext(t)
%     Sext = @(t) sqrt(2*D/dt) * randn();
% 
%     % —— 构造系统右端和雅可比
%     f   = @(t,X) rhs_ext(t, X, a, b, beta, A, omega_k, alpha, mu_s, Sext);
%     j = @(X)  jac(X, a, b, beta, A, omega_k, alpha);  % 不含噪声
% 
%     % —— 计算 LLE
%     lle_vals(k) = LEs(f, j, X0, dt, T, warmup, 20);
% end
% 
% %% —— 绘图 —— 
% figure;
% plot(omega_list, lle_vals, '.-','MarkerSize',6); hold on;
% yline(0,'r--','LineWidth',1.2);
% xlabel('\omega','FontSize',12);
% ylabel('最大李雅普诺夫指数','FontSize',12);
% title(['LLE vs \omega (\alpha = ' num2str(alpha) ', D = ' num2str(D) ')'],'FontSize',14);
% grid on;
% 
% % rhs_ext.m ——— 带外部刺激的忆阻神经元模型（可选 Sext）
% % 用法：
% %   dX = rhs_ext(t, X, a, b, beta, A, omega, alpha, mu_s)
% %   dX = rhs_ext(t, X, a, b, beta, A, omega, alpha, mu_s, Sext)
% %
% % 如果不传入 Sext，则默认 Sext = 0。
% function dX = rhs_ext(t, X, a, b, beta, A, omega, alpha, mu_s, Sext)
%     % 处理可选的 Sext
%     if nargin < 10 || isempty(Sext)
%         s = 0;
%     else
%         % 如果 Sext 是函数句柄，则调用它，否则按常量处理
%         if isa(Sext,'function_handle')
%             s = Sext(t);
%         else
%             s = Sext;
%         end
%     end
% 
%     x = X(1);
%     y = X(2);
%     z = X(3);
% 
%     % 带 Sext 的方程
%     dx = -alpha*((y - x) - (y - x)^3/3 + (a + 3*b*z^2)*x) + s;
%     dy =      (y - x) - (y - x)^3/3 - beta*y + mu_s + A*cos(omega*t) + s;
%     dz =      x + s;
% 
%     dX = [dx; dy; dz];
% end

